import numpy as np
import matplotlib.pyplot as plt

# # # 设置随机种子以确保结果可重复
# # np.random.seed(42)

# # 在[-1, 1]范围内随机生成30000个点
# points = np.random.uniform(-1, 1, (30000, 3))

# # 筛选满足不同条件的点
# mask_x = (points[:, 0] > points[:, 1]) & (points[:, 0] > points[:, 2])
# mask_y = (points[:, 1] > points[:, 0]) & (points[:, 1] > points[:, 2])
# mask_z = (points[:, 2] > points[:, 0]) & (points[:, 2] > points[:, 1])

# points_x = points[mask_x]
# points_y = points[mask_y]
# points_z = points[mask_z]

# # 创建3D图形
# fig = plt.figure(figsize=(12, 10))
# ax = fig.add_subplot(111, projection='3d')

# # 绘制三批不同颜色的点云
# ax.scatter(points_x[:, 0], points_x[:, 1], points_x[:, 2], c='red', alpha=0.6, label='x > y & x > z')
# ax.scatter(points_y[:, 0], points_y[:, 1], points_y[:, 2], c='green', alpha=0.6, label='y > x & y > z')
# ax.scatter(points_z[:, 0], points_z[:, 1], points_z[:, 2], c='blue', alpha=0.6, label='z > x & z > y')

# # 设置坐标轴标签
# ax.set_xlabel('X', fontweight='bold')
# ax.set_ylabel('Y', fontweight='bold')
# ax.set_zlabel('Z', fontweight='bold')

# # 设置坐标轴范围
# ax.set_xlim(-1, 1)
# ax.set_ylim(-1, 1)
# ax.set_zlim(-1, 1)

# # 添加标题
# ax.set_title('Points colored by their maximum coordinate', fontweight='bold')

# # 绘制加粗的彩色坐标轴
# ax.plot([-1, 1], [0, 0], [0, 0], color='red', linewidth=2)    # X轴
# ax.plot([0, 0], [-1, 1], [0, 0], color='green', linewidth=2)  # Y轴
# ax.plot([0, 0], [0, 0], [-1, 1], color='blue', linewidth=2)   # Z轴

# # 在坐标轴正端标注 x、y、z
# ax.text(1.1, 0, 0, "x", color='red', fontsize=12, fontweight='bold')
# ax.text(0, 1.1, 0, "y", color='green', fontsize=12, fontweight='bold')
# ax.text(0, 0, 1.1, "z", color='blue', fontsize=12, fontweight='bold')

# # 移除默认的灰色坐标轴
# ax.xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
# ax.yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
# ax.zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))

# # 调整视角以更好地展示形状
# ax.view_init(elev=20, azim=45)

# # 添加图例
# ax.legend()

# # 显示图形
# plt.show()


# exit()

# def calculate_centroids(points, sample_sizes):
#     return np.array([np.mean(points[:, i, :size, :], axis=1) for i, size in enumerate(sample_sizes)])

# def euclidean_distance(points):
#     return np.sqrt(np.sum(points**2, axis=-1))

# print("Step 1-2: Start setting parameters and initialization")
# # Set parameters
# sample_sizes = np.arange(5, 15)  # Sampling sizes from 5 to 14
# num_experiments = 4000*600
# num_sample_sizes = len(sample_sizes)

# print(f"Parameters set. Sample size range: {sample_sizes}, Number of experiments: {num_experiments}")
# print("Sampling will be done in the range: 0 < x < 1 and 0 < y < 1")

# print("\nStep 3: Generate all random points at once")
# # Generate all random points at once
# max_sample_size = sample_sizes[-1]
# all_points = np.random.random((num_experiments, num_sample_sizes, max_sample_size, 2))
# all_points += np.finfo(float).eps  # Ensure points are strictly greater than 0
# print("Random point generation completed")

# print("\nStep 4-5: Calculate all centroids")
# # Calculate centroids for all experiments and sample sizes at once
# all_centroids = calculate_centroids(all_points, sample_sizes)
# all_centroids = np.transpose(all_centroids, (1, 0, 2))
# print(f"Centroid array shape: {all_centroids.shape}")
# print("First centroid of the first 5 experiments:")
# print(all_centroids[:5, 0])

# print("\nStep 6: Calculate average distances")
# distances = euclidean_distance(all_centroids)
# average_distances = np.mean(distances, axis=0)
# print("Average distance calculation completed")
# print("Average distances for each sample size:")
# for size, dist in zip(sample_sizes, average_distances):
#     print(f"Sample size: {size}, Average distance: {dist:.6f}")

# print("\nStep 7: Plot the curve")
# # Plot the curve
# plt.figure(figsize=(10, 6))
# plt.plot(sample_sizes, average_distances, marker='o')
# plt.title('Average Distance of Centroids from Origin vs. Sample Size')
# plt.xlabel('Sample Size')
# plt.ylabel('Average Distance')
# plt.grid(True)
# plt.show()

# print("Program execution completed")


# def generate_points(shape):
#     # Generate points in the range -1 < x < 1, -1 < y < 1, and x > y
#     points = np.random.random(shape) * 2 - 1  # Generate points in [-1, 1) x [-1, 1)
#     mask = points[..., 0] <= points[..., 1]  # Find points where x <= y
#     points[mask] = points[mask][..., ::-1]  # Swap x and y for these points
#     return points * (1 - np.finfo(float).eps)  # Ensure strict inequalities

# def calculate_centroids(points, sample_sizes):
#     return np.array([np.mean(points[:, i, :size, :], axis=1) for i, size in enumerate(sample_sizes)])

# def euclidean_distance(points):
#     return np.sqrt(np.sum(points**2, axis=-1))

# print("Step 1-2: Start setting parameters and initialization")
# # Set parameters
# sample_sizes = np.arange(5, 15)  # Sampling sizes from 5 to 14
# num_experiments = 10000
# num_sample_sizes = len(sample_sizes)

# print(f"Parameters set. Sample size range: {sample_sizes}, Number of experiments: {num_experiments}")
# print("Sampling will be done in the range: -1 < x < 1, -1 < y < 1, and x > y")

# print("\nStep 3: Generate all random points")
# # Generate all random points
# max_sample_size = sample_sizes[-1]
# all_points = generate_points((num_experiments, num_sample_sizes, max_sample_size, 2))
# print("Random point generation completed")

# print("\nStep 4-5: Calculate all centroids")
# # Calculate centroids for all experiments and sample sizes
# all_centroids = calculate_centroids(all_points, sample_sizes)
# all_centroids = np.transpose(all_centroids, (1, 0, 2))
# print(f"Centroid array shape: {all_centroids.shape}")
# print("First centroid of the first 5 experiments:")
# print(all_centroids[:5, 0])

# print("\nStep 6: Calculate average distances")
# distances = euclidean_distance(all_centroids)
# average_distances = np.mean(distances, axis=0)
# print("Average distance calculation completed")
# print("Average distances for each sample size:")
# for size, dist in zip(sample_sizes, average_distances):
#     print(f"Sample size: {size}, Average distance: {dist:.6f}")

# print("\nStep 7: Plot the curve")
# # Plot the curve
# plt.figure(figsize=(10, 6))
# plt.plot(sample_sizes, average_distances, marker='o')
# plt.title('Average Distance of Centroids from Origin vs. Sample Size')
# plt.xlabel('Sample Size')
# plt.ylabel('Average Distance')
# plt.grid(True)
# plt.show()

# # Bonus: Visualize the sampling distribution
# print("\nBonus: Visualizing the sampling distribution")
# plt.figure(figsize=(8, 8))
# sample_points = generate_points((10000, 2))
# plt.scatter(sample_points[:, 0], sample_points[:, 1], alpha=0.1)
# plt.title('Distribution of Sampled Points')
# plt.xlabel('x')
# plt.ylabel('y')
# plt.axline([0, 0], [1, 1], color='r', linestyle='--', label='x = y')
# plt.legend()
# plt.grid(True)
# plt.show()

# print("Program execution completed")


# def generate_points(shape):
#     # Generate points in the range 0 < x < 1 and 0 < y < 1
#     return np.random.random(shape) * (1 - np.finfo(float).eps)

# def calculate_centroids(points, sample_sizes):
#     return np.array([np.mean(points[:, i, :size, :], axis=1) for i, size in enumerate(sample_sizes)])

# def manhattan_distance(points):
#     return np.sum(np.abs(points), axis=-1)

# def euclidean_distance(points):
#     return np.sqrt(np.sum(points**2, axis=-1))

# print("Step 1-2: Start setting parameters and initialization")
# # Set parameters
# sample_sizes = np.arange(5, 15)  # Sampling sizes from 5 to 14
# num_experiments = 1000000
# num_sample_sizes = len(sample_sizes)

# print(f"Parameters set. Sample size range: {sample_sizes}, Number of experiments: {num_experiments}")
# print("Sampling will be done in the range: 0 < x < 1 and 0 < y < 1")

# print("\nStep 3: Generate all random points")
# # Generate all random points
# max_sample_size = sample_sizes[-1]
# all_points = generate_points((num_experiments, num_sample_sizes, max_sample_size, 2))
# print("Random point generation completed")

# print("\nStep 4-5: Calculate all centroids")
# # Calculate centroids for all experiments and sample sizes
# all_centroids = calculate_centroids(all_points, sample_sizes)
# all_centroids = np.transpose(all_centroids, (1, 0, 2))
# print(f"Centroid array shape: {all_centroids.shape}")
# print("First centroid of the first 5 experiments:")
# print(all_centroids[:5, 0])

# print("\nStep 6: Calculate average Manhattan distances")
# distances = euclidean_distance(all_centroids)
# average_distances = np.mean(distances, axis=0)
# print("Average Manhattan distance calculation completed")
# print("Average Manhattan distances for each sample size:")
# for size, dist in zip(sample_sizes, average_distances):
#     print(f"Sample size: {size}, Average Manhattan distance: {dist:.6f}")

# print("\nStep 7: Plot the curve")
# # Plot the curve
# plt.figure(figsize=(10, 6))
# plt.plot(sample_sizes, average_distances, marker='o')
# plt.title('Average Manhattan Distance of Centroids from Origin vs. Sample Size')
# plt.xlabel('Sample Size')
# plt.ylabel('Average Manhattan Distance')
# plt.grid(True)
# plt.show()

# # Bonus: Visualize the sampling distribution
# print("\nBonus: Visualizing the sampling distribution")
# plt.figure(figsize=(8, 8))
# sample_points = generate_points((10000, 2))
# plt.scatter(sample_points[:, 0], sample_points[:, 1], alpha=0.1)
# plt.title('Distribution of Sampled Points')
# plt.xlabel('x')
# plt.ylabel('y')
# plt.grid(True)
# plt.show()

# print("Program execution completed")


# def generate_points(shape):
#     # Generate points in the range 0 < x < 1 and 0 < y < 1
#     return np.random.random(shape) * (1 - np.finfo(float).eps)

# def calculate_centroid(points):
#     return np.mean(points, axis=0)

# def euclidean_distance(point):
#     return np.sqrt(np.sum(point**2, axis=-1))

# print("Step 1-2: Start setting parameters and initialization")
# # Set parameters
# sample_sizes = np.arange(5, 25)  # Sampling sizes from 5 to 14
# num_experiments = 100000
# num_sample_sizes = len(sample_sizes)

# print(f"Parameters set. Sample size range: {sample_sizes}, Number of experiments: {num_experiments}")
# print("Sampling will be done in the range: 0 < x < 1 and 0 < y < 1")

# print("\nStep 3: Generate all random points and calculate centroids")
# # Generate all random points and calculate centroids
# overall_centroids = np.zeros((num_sample_sizes, 2))
# for i, size in enumerate(sample_sizes):
#     points = generate_points((num_experiments, size, 2))
#     overall_centroids[i] = calculate_centroid(points.reshape(-1, 2))

# print(f"Overall centroid shape: {overall_centroids.shape}")
# print("Overall centroids:")
# print(overall_centroids)

# print("\nStep 4: Calculate Euclidean distances")
# distances = euclidean_distance(overall_centroids)
# print("Euclidean distance calculation completed")
# print("Euclidean distances for each sample size:")
# for size, dist in zip(sample_sizes, distances):
#     print(f"Sample size: {size}, Euclidean distance: {dist:.6f}")

# print("\nStep 5: Plot the curve")
# # Plot the curve
# plt.figure(figsize=(20, 6))
# plt.plot(sample_sizes, distances, marker='o')
# plt.title('Euclidean Distance of Overall Centroid from Origin vs. Sample Size')
# plt.xlabel('Sample Size')
# plt.ylabel('Euclidean Distance')
# plt.grid(True)
# plt.show()

# # Bonus: Visualize the sampling distribution and overall centroids
# print("\nBonus: Visualizing the sampling distribution and overall centroids")
# plt.figure(figsize=(8, 8))
# sample_points = generate_points((10000, 2))
# plt.scatter(sample_points[:, 0], sample_points[:, 1], alpha=0.1, label='Sample Points')
# plt.scatter(overall_centroids[:, 0], overall_centroids[:, 1], color='red', s=100, label='Overall Centroids')
# plt.title('Distribution of Sampled Points and Overall Centroids')
# plt.xlabel('x')
# plt.ylabel('y')
# plt.legend()
# plt.grid(True)
# plt.show()

# print("Program execution completed")

# def generate_points(shape):
#     # Generate points in the range 0 < x < 1 and 0 < y < 1
#     return np.random.random(shape) * (1 - np.finfo(float).eps)

# def calculate_centroid(points):
#     return np.mean(points, axis=0)

# print("Step 1-2: Start setting parameters and initialization")
# # Set parameters
# sample_sizes = np.arange(5, 15)  # Sampling sizes from 5 to 14
# num_experiments = 10000
# num_sample_sizes = len(sample_sizes)

# print(f"Parameters set. Sample size range: {sample_sizes}, Number of experiments: {num_experiments}")
# print("Sampling will be done in the range: 0 < x < 1 and 0 < y < 1")

# print("\nStep 3: Generate all random points and calculate centroids")
# # Generate all random points and calculate centroids
# all_centroids = []
# for size in sample_sizes:
#     points = generate_points((num_experiments, size, 2))
#     centroids = np.array([calculate_centroid(p) for p in points])
#     all_centroids.append(centroids)

# print("Centroid calculation completed")

# print("\nStep 4: Plot all centroids")
# plt.figure(figsize=(12, 10))

# # Use a color map to assign different colors to different sample sizes
# colors = plt.cm.rainbow(np.linspace(0, 1, num_sample_sizes))

# for i, (size, centroids) in enumerate(zip(sample_sizes, all_centroids)):
#     plt.scatter(centroids[:, 0], centroids[:, 1], color=colors[i], alpha=0.1, 
#                 label=f'Sample size: {size}')

# plt.title('Distribution of Centroids for Different Sample Sizes')
# plt.xlabel('x')
# plt.ylabel('y')
# plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
# plt.grid(True)
# plt.tight_layout()
# plt.show()

# # Bonus: Visualize the sampling distribution
# print("\nBonus: Visualizing the sampling distribution")
# plt.figure(figsize=(8, 8))
# sample_points = generate_points((10000, 2))
# plt.scatter(sample_points[:, 0], sample_points[:, 1], alpha=0.1, label='Sample Points')
# plt.title('Distribution of Sampled Points')
# plt.xlabel('x')
# plt.ylabel('y')
# plt.legend()
# plt.grid(True)
# plt.show()

# print("Program execution completed")